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The reader can find in the literature a lot of different techniques to study the dynamics of a given 
system and al so, many suitable num erical integrators to compute them. Notwithstanding the 
recent work of Mafhone et a/.l . l2011al | for mappings, a detailed comparison among the widespread 
indicators of chaos in a general system is still lacking. Such a comparison could lead to select the 
most efficient algorithms given a certain dynamical problem. Furthermore, in order to choose the 
appropriate numerical integrators to compute them, more comparative studies among numerical 
integrators are also needed. 



This work deals wit h both problem s . We first extend the work of (Maffione et Idl . l2011aj fo 



mappings to the 2D |Henon fc Heilesl . ll964l| potential, and compare several variational indicators 
of chaos: the Lyapunov Indicator (LI); the Mean Exponential Growth Factor of Nearby Orbits 
(MEGNO); the Smaller Alignment Index (SALI) and its generalized version, the Generalized 
Alignment Index (GALI); the Fast Lyapunov Indicator (FLI) and its variant, the Orthogonal 
Fast Lyapunov Indicator (OFLI); the Spectral Distance {D) and the Dynamical Spectras of 
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Stretching Numbers (SSNs). We also include in the record the Relative Lyapunov Indicator 
(RLI), which is not a variational indicator as the others. Then, we test a numerical technique 
to i ntegrate Ordinary D ifferential Equations (ODEs) based on the Taylor method implemented 
by |jorba fc Zoul . l2005j (called t aylor), and we compare i ts performance with other two well- 



known efficient integrators: the (Prince fc Dormandl . 1981 implementation of a Runge-Kutta 
of order 7-8 (D0PRI8) and a Bulirsch-Stoer implementation. These tests are run under two 
very different systems from the complexity of their equations point of view: a triaxial galactic 
potential model and a perturbed 3D quartic oscillator. 

We first show that a combination of the FLI/OFLI, the MEGNO and the GALl2Ar succeeds in 
describing in detail most of the dynamical characteristics of a general Hamiltonian system. In 
the second part we show that the precision of taylor is better than that of the other integrators 
tested, but it is not well suited to integrate systems of equations which include the variational 
ones, like in the computing of almost all the preceeding indicators of chaos. The result which 
induces us to draw this conclusion is that the computing times spent by taylor are far greater 
than the times insumed by the D0PRI8 and the Bulirsch-Stoer integrators in such cases. On 
the other hand, the package is very efficient when we only need to integrate the equations of 
motion (both in precision and speed), for instance to study the chaotic diffusion. We also notice 
that taylor attains a greater precision on the coordinates than either the D0PRI8 or the 
Bulirsch-Stoer. 

Keywords: Chaos Indicators - Numerical Integrators - Variational Equations - ODEs - Hamil- 
tonian Systems 



1. Introduction 



In order to study a Hamiltonian system, it is essential to count on both fast and easy-to-compute techniques 
to determine the chaotic or regular nature of its orbits. 

For two degrees of freedom (hereafter d.o.f.), a technique based on a graphical treatment, like the 
Poincare Surfaces of Section, should s uffice to understa n d such a low dim ensional problem. Though gen- 
eralized to three dimensional systems ( [Froeschlel . 197Clf |. Froeschle . 1972f |). the graphical treatment shows 
fundamental constraints. Therefore, non graphical methods seem to be the logical alternative to deal with 
systems of higher dimensionality. 

The concept of exponential divergence is thus introduced, and it is one of the two wides pread ideas i n 
chaos detection. The introduction of the L yapunov Characterist ic Exponents (LCEs; see e.g. [ Reichll . 1992l |) 



as well as its numerical implementation ( jBenettin" were a breakthrough contribution to the 

field. Notice that the time of integration is limited, so we are a ble to reach just trunc a ted approximati ons 
of the theoretical LCE s, i.e. Lyapunov In dicators, hereafter LI ( [Benettin et al. . 1976], Froeschle . 1984 ] or 



Tancredi et all l200l[ | and [Skokosl . l201Clf | for reviews on the subject). Therefore, a drawback in the com- 
putation of the truncated values of the LCEs is their very slow convergence. Nevertheless, a large number 
of variational techniques i mprove the Li's perform ance, such as: the Mean Exponential Growth factor of 
Nearby Orbits (MEGNO; [Cincotta fc SimdI . l2nonfh : the Smaller Ahg nment Index ( S ALL jSkokosI, [joOlllgl) 



and its generalized vers ion, the General ized Alignment I ndex (GALL [Skokos et al\ . 120071 ) I) ; the Fast Lya- 



punov Indicator (FLI; (Froeschle eTaH [l99 7]. | Froeschl e et al\ . Il997al |^ and its vari ant, the Orth og;onal 



^ ^_ J ^ ■&MMMin fcM • J ' L — , I J—I/ -> ^ . 

Fast Lyapunov Indicator (OFLI; [Fouchard et aUl2002l ]): the Spectral Distance {D; [Voglis et all 119991 ]) 



200; 



^ See also Cincotta et al\. 
iLemaitre e t Hinse et al 



20041. 



Gozdziewski et all 

¥r 



Skokos et adbOOSl]. [Manos fc Athanassoulal. I20TI 
Froesc hle fc Lega'.'l998l. 'Froeschle fc Lega'. '2000*1. ^ 
|2006|] . [Paleari et a/.. ,2008,] . Todorovic et a/.. . .2008,] . ^Lcga et al 



l2005l |. [Gavon fc Boig . |2008| ]. 



Giordan o fc Cincottal. _ .. . _ 

2010], MafEone et aL','201ll, 'Compere et a/.','201l|] 

Skokos et al, 2004], Szell et al, 2004], Bountis fc Skokos. 2006], Carpinterol, I2OO8I ] . [Antonopoulos et aZ.I. [2OI0I ] . 



See also 

^See also ^^^^^^^^^^^^ 

"See also [Froesc hle fc Lega'.'l998^. 'Froeschle fc Lega'.'2000']. ^Lega fc Froeschldl200lj . [Cuzzo et fflZl.l2002t l. [Froeschle fc Legal . 



.2010|] . 



Comparative study of variational chaos indicators and ODEs ' numerical integrators 3 



and the Dynamical Spectras of Stretching Nu mbers (SSNs: jVoglis &: Contopoulos . 1994 1). Finally, we 
include the Relative Lyapunov Indicator (RLI; (Sandor et al. . 20001 ^ . which is not based on the evolution 
of the solution of the first variational equations as the others, but on the evolution of two different but 
very close orbits. 

These are just a sample of all the variational Chaos Indicators (CIs) widely used nowadays. They are 
the ones selected to be used in the forthcoming investig ation. However, the record could indeed continue 



20081); th e Orthogo- 



to include the Average Power Law Exponent (APLE: iLukes-Gerakopoulos et a/.l. . „ . . „ 

nal F ast Lyapunov Indicator of order two (0FLl2^; [Barrio . 2005], Barrio et a/.l . 12009.] . [Barrio et al. 
201II which is a secorid ord er variational method; the Dy n amical Spectras of Helic it y or Twist Angles 
rIContopoulos fc Voglisl . ll996l ]. [Contopoulos fc Voelisl. Il997ll. iK-oeschle fc Legal . Il998l |. [Lega fc Froeschlel 
1998 1]). Moreover, evolutionarv algorithms ( Petalas et al. . 20091 ]) are also interesting alternatives. 



The other widespread idea in the field of CIs is the analysis of some particular qua ntities (e.g. the 
frequency) of a single orbit. In this group we find the Frequency M ap Analysis (FMA: (Laskarl. [l 990l^) 
and its variant, the Frequency Modified Fourier Transform ( FMFT; ISidlichovskv fc: Nesvorni . 19961]) o r 
other alternatives using F a st Fourier Transforni tech niques ( |Michtchenko. T. fc Ferraz-Mello. S. . 1995l |. 



(Ferraz-Mello et ali^m^ . (Michtchenko et a?.! . [2OI0I ]). among others. Nevertheless, we leave such family 



of indicators for a future comparison. 

In order to compute the above variational CIs efficiently, it is fundamental to appropriately recognize 
numerical techniques for the integration of Ordinary Differential Equations (ODEs hereafter), which is the 
goal of the second part of this review. That is, in most of the problems that arise from the investigation of 
dynamical systems, as well as in every realistic model, the corresponding differential equations have to be 
solved by means of numerical techniques. Its computational implementation should be as fast and precise 
as possible, in order to minimize both the numerical error and the computing time. 

Among all numerical integrators, there are those which use recurrent power series in order to compute 
the solutions wit h a considerab lv high precis i on, an d if we take advantage of the recurrence, the computing 
time is reduced 
Taylo r meth o d ( 



Sitarskil. [19791 fc iSitl^. Il989l]). One of t he mos t used p ower series integrators is the 



MontenbrudIll99ll]. |Montenbrud3. Il992| . iGofenl. 1199^1. |Gozdziewski fc Maciejewski 



199.51 ]. (SimdI . l200()l ]. [Simol . 120081 ] fc [Gerlach fc Skokosl . l201ol ]V whose usefulness is at all means clear 



particu lar implement ation of which will be tested in this work. 

In [jorba &: Zou . 2005 ] the authors presented an integration package (coded in C) based on the Taylor 
method for numerical integration of ODEs. They showed that this package turned out to be a very efficient 
integrator (both in speed and precision) when applied to the Lorentz system, the forced pendulum model 
and the Restricted Three-Body Problem. 

The reason for selecting this particular package for our numerical experiments is that it implements 
several features which makes it a very versatile and nearly automatic package. This is due to the fact that 
the software encompasses the necessary subro utines to obtain a co mple t e inte g rator, such as t he iiiclusion 



of the so-called automatic d ifferentiation (iBeda et al. . 1959.]. iWengert . 1964 1. Moore . 1966ll . Rail . 1981 ] 



Griewank fc Corhsd . [l99ll ]. [Bischof et a"n . ll992l ]. [Berz et a/.l . 1l996l ] and [Criewankl . l200oj )^ an algorithm 



to choose the order of the power series p and the step size h, the possibility to incorporate customized 
mathematical routines that use higher precisions than the standard double precision, and even a FORTRAN 
wrapper to call the integrator from a main program coded in FORTRAN 77. 

Taking into account those advantageous features of Jorba & Zou package, we have decided to test their 
implementation of Taylor method by analizing its performance when applied to three different dynamical 
problems for which we have compared its efficiency -in both computing time and precision- against that 
of other two widely used integrators. 

This work comprises two parts: in the first part of the paper, we make a comparative study of variational 
CIs including the LI, the RLI, the MEGNO, the D and the SSNs, the FLI and the OFLI, the SALI and 
the GALIs. In Section [21 we briefly introduce the above-mentioned CIs and in Section [3l we perform some 



^See also 
^See also 
See also 



Contopoulos fc yogliEl,[l99^. [Contopoulos fc Voglislll997ll [Contopoulos et aZ.l[l997at . [Voglis et fflZ.l[l996 

Szell et aZ.l [20041. [Sandor et al)}2m'^. 'Sandor et a/.l'2007|. 

Laskar et a/.L [199^ 1 . [Papaphilippou fc Laskar . 1996.] . ^Papaphilippou fc LaskaJ . [l998| |. 
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experiments on the simple model of Henon &: Heiled . Il964l | (hereafter HH). There, we deal with the CI's 
final values and time evolution curves. In Section 13.1^ we study the reliability of the CI's thresholds in 
order to separate chaotic from regular motion, and also the level of dynamical information provided by the 
CI's final values. Then, in Section [3.21 we examine the capability of identifying periodicity and instabilities 
by means of the CI's time evolution, and the sensitivity of the CIs on initial deviation vectors (hereafter, 
i.d.v.). Finally, we discuss the computing times and the advantage of the saturation times ( see Sectionl3.3l) . 
Our p urpose in this first part of the reported research is to extend the results shown in (Maffione et al\ . 



2011al ] to a Hamiltonian flow for mappings. 



Then, in the second part of the paper, the comparison among numerical techniques for the resolution 
of ODEs in order to compute the CIs just mentioned is presented. Some of the features of Jorba & Zou 
package will be briefly described in Section HI We introduce the results of testing the latter package over 
three different dynamical problems as follows: the first two problems (Section [5] and Section [6]) concern 
the computation of two dynamical indicators, i.e. the MEGNO (Section 15.11 and Section 16. ip and the LI 
(Section 15.21 and Section l6.2p . Both problems require the rather accurate integration of the equations of 
motion along with their first variationals and, in the first case, some additional differential equations must 
be included. The third problem under consideration is the solely integration of the equations of motion 
(see Section [5.31 and Section [6. 3p . Though this is a rather simple problem, it has many useful applications, 
such as the study of chaotic diffusion in phase space, where we need to follow accurately the coordinates 
on phase space over rather long time intervals. For the above-mentioned problems, we have chosen two 
particular Hamil tonian systems: the potential corresponding to an auto-consistent model of a triaxial 
elliptical galaxy Muzzio et al . 2005l |. which actually involves complicated expressions for the ODEs, and 



a rather simpler one, a perturbed three dimensional quartic oscillator. In Section [7] we only compute the 
equations of motion for the three integrators to test the precision of the coordinates of phase space, and 
compare which one preserves better the position on the orbit. Finally, we discuss the results in Section [51 



2. Selected CIs revisited 

The aim of the present section is to introduce the CIs that we will use in the forthcoming investigation. 
However, for the sake of simplicity, we describe only the main properties defined for flows. 



2.1. The LI and the MEGNO 

Consider a continuous dynamical system defined on a differentiable manifold S, where $*(x) characterizes 
the state of the system at time t and being the state of the system at time t = 0, x(0). Therefore, the state of 
the system after two consecutive time steps t and t' will be given by the composition law: = o . 

Moreover, the tangent space of x maps onto the tangent space of $*(x) according to the operator 
dx^* and following the rule w(t) = (ix^*(w(0)) where w(0) is an i.d.v. The action of such operator at 
consecutive time intervals satisfies the equation: 

If we suppose that our manifold S has some norm denoted by || • ||, we can define the useful quantity: 



||dx**w|| 
||w| 

called "growth factor" in the direction of w. 

Let us now say that H{p, q) with p, q G be an A^-dimensional Hamiltonian, that we suppose 
autonomous just for the sake of simplicity. Introducing the following notation: 

X = (p,q) G M^^, f(x) = i-dH/dq, dH/dp) G M^^, 
the equations of motion can be written in a simple way like 
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x = f(x). (1) 

Let 7(xo; t) be an arc of an orbit of the flow ([T|) over a compact energy surface: Mh C M^^, = {x : 
ff(p, q) = h} with h = constant, then 



7(xo;t) = {x(t';xo) : xq € A^, < < t}. 



We can gain fundamental inf ormation abo ut the Hamiltonian flow in the neighborhood of any orbit 7 
through the largest LCE (ILCE, [Skokosl . [2ninl |) defined as: 



X[7(xo;0] = lim -lnAt[7(xo;t)] 

t— >-oo t 



(2) 



with 



A,[7(xo;f)] = " I „ " 
||w|| 

where ||(i^$*w|| is an "infinitesimal displacement" from 7 at time t. The fact that the ILCE (and its 
truncated value, the so-called LI= limt^j^ ]■ In At[7(xo; t)] for T, finite) measures the mean exponential 
rate of divergence of nearby orbits is clearly understood when written Eq. Q in an integral fashion: 



1 /■* 

x[7(xo;t)] = hm - / 

t^oo t Jq 



\d-./^^'w\\ 



(3) 



where the bar denotes time average. 

Now we are in condition to introduce the MEGNO f [Cincotta fc Simol . [200ol | ) . y[7(xo;t)], through the 
expression: 



t Jq w|| 

which is related to the integral in Eq. ([3]); i.e., in case of an exponential increase of ||(i^$*w||, ||(i^$*w|| = 
||w|| • exp(xt), the quantity y[7(xo;i)] can be considered as a weighted variant of the integral in Eq. ([3]). 
Listead of using the instantaneous rate of increase, Xi we average the logarithm of the growth factor, 
In(At) = xt- 

Introducing the time average: 



y[7,(xo;t)] = - / y[7,(xo;t')]dt', 



(5) 



we notice that the time evolution of the MEGNO can be briefly described in a suitable and unique 
expression for all kinds of motion. Indeed, its asymptotic behavior can be summarized in the following 



way: y[7(xo;t) 



a^t + b-y, where 



X7/2 and bj ^ for irregular, stochastic motion, while = 



and 6^ ~ 2 for quasi-periodic motion. Deviations from the value 6^ ~ 2 indicate that 7 is close to some 
particular objects in phase space, being 6^ < 2 or 6^ > 2 for s table periodic orbits (re sonant elliptic tori), 
or unstable periodic orbits (hyperbolic tori), respectively (see Cincotta &: Simd . 
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2.2. The FLI and the OF LI 

The FLI (see [Froeschle etdl . Il997l ]. (Froeschle et ali Il997a 
timately related to the ILCE (and the MEGNO 



see 



Froeschle h Legal . l2006l | ) is a quantity in- 



Mestre et ali 1201 ll ]). which is able to distinguish 



b etween chaotic (w e akly ch aotic) and regular inotion and also, between resonant and non resonant motion 
( [Froeschle fc Legal . [2nnnl |. [Lega fc Froeschlel . l2nni[ |) using only the first part of the computation of the 
ILCE. 

On an A^-dimensional flow, we follow the time evolution of N unit i.d.v. Therefore, the FLI at time t 
is defined by the highest norm among the unit i.d.v. of the basis, as follows: 



FLI{t) = supt [||w(t)i||, ||w(t)2||, . . . , ||w(t)jv||] . 

For further details, refer to jFroeschle et al , 1997al ]. 

For both kinds of motion, the FLI tends to infinity as the time increases, with completely different 
rates, though (it behaves linearly for regular moti o n and grows exponentially fast for chaotic motion). 

In the case of the OFLI (see (Fouchard et al , 2002]), we take the component orthogonal to the flow 
for each unit i.d.v. of the basis at every time step: 



OFLI{t) = supt w{t)j;,w{t) 



)2 ' ■ 



,w{t) 



N 



This modification makes the OFLI a CI that can easily distinguish periodicity among the regular 
component. The OFLI for periodic orbits oscillates around a constant value, while for quasiperiodic motion 
and chaotic motion has the same behavior as the FLI. 



2.3. The SSNs and the D 

The stretching number Sj is defined as (s ee Contopoulos &: Voglis . 1996], [Contopoulos et al. . 1997a | 
Contopoulos fc Voghsl . ll997l |. (Voglis et all 11997 * 



]dx* 



t+i-5t 



(w(0))] 



(^t^"" ]dx**+(*-i)-^*(w(0))l 

where dx^*'^*''^'(w(0)) = w(t + i ■ 6t), i.e. the tangent vector at time t + i ■ 6t, with i = 1,2, . . . ,p and p is 
the number of steps of integration 6t (in case of mappings, the 6t = 1) in which the interval of integration 
is divided. Therefore, the serie of Si with i oo converges to the ILCN. 

Another way to obtain dynamical information is analyzing the spectra of the Si, i.e. the SSNs. However, 
when the sample to study is very lar ge, the single analysis of the orbits by means of the SSNs is no longer 
reliable. Thus, in IVoglis et all Il999t ] the authors introduce the D. 

The SSNs are defined by the probability density of the values s the Si can take: 



S{s) 



dZ{s, s + ds) 
Z ■ ds 



with Z the total number of si and dZ{s, s + ds) the number of Si in the interval {s, s + ds). 
There are certain properties of the SSNs histograms that are worth mentioning: 

(1) The dynamical spectra S{s) is invariant regarding the initial condition (hereafter, i.e.) along the same 
orbit (invariant with respect to time). 

(2) The dynamical spectra S{s) is invariant regarding the i.e. in a connected chaotic domain (invariant 
with respect to space). 

Finally, we define the D like: 



D^ = Y,[Si{s)-S2{s)f-ds, 

s 
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where Sj{s) is the normahzed number of Si corresponding to the unit i.d.v. Wj(0) which have values in the 
interval (s, s + ds) (it is necessary to follow the evolution of two unit i.d.v. for each i.e. to compute the D). 

Therefore, based on the properties of the SSNs, it is straighforward to see that if the orbit is regular, 
the D tends to a constant non-zero value. If the orbit is chaotic, the D decreases approaching to zero. 

2.4. The SALI and the GALI 



In [Skokosl . l200l[ | the author introduces the SALI defining the alignment indices. The parallel alignment 



index: 

d_ = II wi — W2II 

and the anti-parallel alignment index: 

d+ = ||wi + W2I 

where the || • || is the usual euclidean norm. If the deviation vectors are normalized, then in case of chaotic 
motion: (i_ — )• and (i+ ^ 2 (the deviation vectors eventually become aligned with the same direction) or 
d- — )• 2 and (i+ — ?> (with opposite directions). If the motion is regular, d- and d+ oscillate within the 
interval (0,2). Therefore, the SALI is defined as the smaller of those indices: 

SALI{t) = min(d-^-,d^). 
The behavior of the SALI for chaotic orbits follows an exponential law: 



SALI{t) oc e 



-{Xi-X2)-t 



shown in [Skokos et al\ . l2004i |. where Xi with i = 1,2 are the two largest LCEs of the orbit. If the orbit is 



regular , the SALI os c illate s within the interval (0,2). 



In [Skokos et all l2007l | the authors generalize the SALI introducing an alternative way to compute it. 



They evaluate the quantity: 

P{t) = d+-d-, 

for each time step, where 

Pit) „ 

= ||Wi A Vi^2|| 

with wi A W2 the wedge product of two deviation vectors. This represents the area of the parallelogram 
formed by the two deviation vectors wi and W2. In fact, the wedge product can be generalized to represent 
the volume of a parallelepiped formed by the deviation vectors wi, W2, . . . ,Wfc with 2 < k < 2N, of an 
A^-d.o.f. Hamiltonian system, or a 2A^-dimensional symplectic mapping. Therefore, the GALI^, is defined 
as the volume of the /c-parallelepiped formed by the k initially linearly independent unit deviation vectors 
Wi{t) = ||^'|*||| , i = 1,2, . . . , k and can be computed through the wedge product as follows: 

GALh{t) = \\wi{t) A W2{t) A ... A Wk{t)\\. (6) 

From Eq. ([6]) it is easy to see that if GALIfc(t) tends to zero, this implies that the volume of the paral- 
lelepiped having the unit i.d.v. as edges also shrinks to zero, as at least one of the deviation vectors becomes 
linearly dependent on the remaining ones. Thus, the behavior of the GALIjt(t) follows the exponential law: 



GALIkit) (X e 



-Kxi-X2)+{xi~xs)+---+{xi-Xk)]-t 
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where Xi; X2, • • • , Xfc are the k largest LCEs of the orbit. 

On the other hand, if GALIfc(t) is different from zero as t increases, the hnear independence of the unit 
i.d.v. involved becomes clear. This second case corresponds to regular motion. Regular orbits move on tori 
of dimension M < A^. In general, the dimension of the tori is because there are no resonances involved. 
The number of resonances diminishes the dimensionality of the manifold and, for periodic motion, the orbit 
mov es on an invariant curve, i . e. a in a nifold of dimension o ne: M = 1. The beh av ior of the GALI^. is given 
by f |Christodoulidi fc Bountisl . l2006l | . |Skokos et adl2007l ]. (Skokos et adl2008l ]. (Manos et adl201ll ]): 



constant if 2 < k < M 
GALIk{t)oc{ ifM<k<2N-M 

if2N -M <k<2N. 



2.5. The RLI 

The definition of the RLI is straightforward (see Sandor et al. . 2004l |). Consider the LI for a given i.e. xq 
and after j steps of integration: L/(xo; j). The RLI is thus defined according to the LI difference for the 
"base" orbit and its "shadow" : 



AL/(xo;i) = ||L/(xo + Ax; j) - L/(xo; j)|| 

, where xq and xq + Ax are the i.e. for the base orbit and its shadow, respectively. 

The initial separation between the base orbit and its shadow, i.e. ||Ax||, is a free parameter. 

It is advisable to smooth the evolution curve of the LI for the base orbit to eliminate fast fluctuations. 
Thus, we define the RLI as the total average each step of the integration: 

, t/St 

RLI{t) = (AL/(xo))t = -5^AL/(xo,z-5t), 

i=l 

5t being the step of integration. 

RLI values for chaotic motion are several orders of magnitude higher than those associated with regular 
motion. 



3. A CIs' function for a Hamiltonian flow 

There are important characteristics that make a CI suitable for a given study. In order to select the adequate 
CIs, it is useful to compare their performances according to these main aspects. Therefore, it is essential 
to study the resolving power, i.e. the CPs capability to distinguish chaotic orbits from regular orbits, 
to describe the levels of chaoticity and the levels of regularity and to identify periodicity and stickiness. 
Moreover, we need to analyze the integration time that the CI requires to identify clearly chaotic and 
regular orbits, i.e. the so-called speed of convergence. Finally, it is also important to study the CPU or 
computing times of the CIs, as we will see in Section 13.31 

The resolving power and the speed of convergence tell us about the possibility of gathering detailed 
dynamical information in short integration times. The CPU times tell us about the limitations of the 
numerical implementation of the algorithms to compute the CIs. 

There are two ways of collecting this dynamical information. The first one is using the time evolution 
of the CI. The second one is by means of the CI's final value or its value at the end of the integration 
process. The information provided by the CI's time evolution is much more detailed than only a record 
of the final value. However, the study of large samples of CI's time evolution curves is not as efficient as 
the study of the corresponding final values despite the difference in the amount of information provided by 
each method. Thus, the use of the CI's time evolution is recommended to analyze small samples of orbits 
or to study particular interesting cases. On the other hand, we use the CI's final values to study extended 
regions of the phase space with large samples of orbits. 
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The aim of the work done by Maffione et aLl . l2011al | was finding a "CIs' function" (hereafter, CIsF) for 
mappings, which means a small package of CIs that represents the most efficient way to gather dynamical 
information. There, the authors dealt with many CIs' characteristics previously mentioned, such as the 
identification of chaotic, regular and sticky orbits, the sensitivity to hyperbolicity and stability levels, the 
speed of convergence and the reliability on the thresholds. However, they did not consider their capability 
to identify periodic orbits, their performances' dependence on i.d.v. and their computing times. In this 
section, we do not only extend their work considering a Hamiltonian fiow, the 2-d.o.f. HH potential, but 
we also increase the number of variational CIs in the comparison, including the OFLI and the GALI^ with 
2 < < 4. 

In Section fS.l.ll we consider the robustness and reliability of the CIs' thresholds in distinguishing chaotic 
from regular motion. As we are interested in studying extended regions of the phase space, following the 
time evolution of the CIs for every single orbit proves unsuitable. Therefore, it is very important to count on 
thresholds that make a confident distinction between chaotic and regular motion. Final values are required 
for big samples of orbits whereas the time evolution of the CIs is the most efficient way to study small 
samples of orbits or individual orbits. In Section 13.1.21 we deal with the CIs' final values and in Section 
3.21 we deal with the CIs' time evolution. We analyze the capability of the CIs to identify periodicity 
(Section I3.2.ip and the sensitivity of the methods on dynamic instabilities (Section I3.2.3P using the CIs' 
time evolution, and the CIs' dependence on the i.d.v. (Section I3.2.2p . Finally, in Section 13.31 we briefiy 
discuss the computing times of the different CIs. 

All the computations in the investigation were done using the following configuration. Hardware: CPU, 
2 X Dual XEON 5450, Dual Core 3.00GHz; M.B., Intel S5000VSA; RAM, 4GB(4xlGB), Kingston DDR-2, 
667MHz, Dual Channel. Software: gfortran 4.2.3. 

The computation of the CIs was accomplished by means of a single code we are developing, the so-called 
La Plata Variational Indicators Code (LP-VIcode). The code uses the D0PRI8 routine. 

We use the following configuration for the rest of Section [3] for the computation of the experiments 
unless stated otherwise: we record the data every 10'^ units of time (hereafter u.t.) u ntil a final integratio n 
time of 10^ u.t. The initial separation taken for the calculation of the RLI is 10" r ISandor eA. a,l\ . \2m^ ^. 
The D is computed at intervals of 10^ time steps, and the number of cells considered for the generation 
of the histograms for the SSNs is 10^. The basis of unit i.d.v. (or simply i.d.v.) is the canonical basis. Let 
us recall that it is important to keep the same i.d.v. for the wh ole sample along the sin gle experiments 
because some dependency of the CIs on the i.d.v. may be found ( [Proeschle &: Lega . [2003]). In Section I 
we deal with such dependence. 



3.1. The CIs' final values 

We start considering the reliability on the thresholds because they are essential to study the resolving 
power and the speed of convergence by means of the CIs' final values in Section 13.1.21 

3.1.1. The reliability on the thresholds 

The following study is undertaken on the 2D HH potential on an energy surface defined by /i = 0.118. 
We adopt a sample of 1.25751 x 10^ i.e. taken in the region defined by a; = 0, y € [—0.1 : 0.1] and 
Py G [-0.05 : 0.05]. 

The well-known 2D HH system is described by the Hamiltonian: 

where x,y,px,Py are the usual phase space variables. 

We applied the time-dependent threshold of the LI (Tabled]) with t as the time. It is known that an 
empirical adjustment of the Li's theoretical threshold is strongly advisable, in particular when studying 
large samples of orbits. Yet, our choice here is to avoid this fine tunning and to consider the raw theoretical 
estimation for the sake of a fair comparison. The critical value used for the RLI (Table [l]) was computed 
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following (Sandor et al . 2004 ] and the remarks discussed in [Maffione et al . 2011a]. The way to determine 
the time-dependent threshold for the D is not known yet, thus it is not considered in this section. 

In the case of the MEGNO, the threshold is a fixed value (Table [1]) which also needs an empirical 
adjustment. However, as we did for the LI, we use the theoretical fixed value. 

For the SALI there are two thresholds commonly used, namely, 10~^^ and 10~^ (see e.g. Skokos et al\ 



20041 ] In between, the orbits are called "sticky chaotic". Nevertheless, we consider them also as chaotic 
orbits. Then, the threshold to be analyzed is lO"'^, which separates regular orbits from chaotic and sticky 
chaotic ones. Once again, this is done to avoid any advantages in taking more than one critical value. 

The threshold a ssociated with the FLI or the OFLI (sometimes used with two thresholds also, see 
(Paleari et all 2008*]) is time-dependent and follows the formulae presented in Tabled) 

Finally, we define the thresholds for the GALIs. The HH potential is a 2-d.o.f. Hamiltonian system 
and we can compute three GALIs: the GALI2, the GALI3, and the GALI4. In order to establish their 
thresholds for large samples of orbits, we considered the formulaes for chaotic m otion first: GALI2 oc e~^^* , 
GAL Is PC e~^^^*, and GAL I4 oc e~^^^*, where yi is the ILCE of the orbit (see Christodouhdi Bountid . 



20061 ] ■ [Skokos et all . l2007l ] and [Manos et all . I2OIII ] or Section [23] of this paper for further details). The 



ILCE was approximated by the LI, which has the threshold ln{t)/t (see Tabled]). Then, the thresholds 
for the GALIs to study large samples of orbits are given in Table [TJ We are aware that the behaviors of 
the GALIs for regular motions change with the i.d.v. initially tangent to the torus. However, the selected 
time-dependent thresholds remain a good approximation to separate regular from chaotic motion. 



CI 


Threshold 


LI 


ln{t)/t 


RLI 


10-^" 


MEGNO 


2 


SALI 


10-^ 


FLI/OFLI 


t 


GALI2 


t-' 


GALI3 


t-'' 


GALI4 


t-^ 



In Maffione et all ]2011al ] the authors tested the reliability of the thresholds of the CIs by evaluating 
the percentage variation of the chaotic component according to a change in their thresholds by ±1%. This 
change emulates a fine tunning of the critical value. Here, we repeat the experiment in the HH potential 
with a slight modification and we find similar results to those showed in the above mentioned work for 
mappings. 

First, we recorded the time evolution of the difference between the percentages of chaotic orbits ac- 
cording to a change in the threshold by ±1% for the LI, and used it to measure the robustness of the 
other CIs' thresholds. That is, we computed the ratios between the CIs and the Li's percentage variation 
of the chaotic component according to the change above-mentioned. If the ratios are above one, it means 
that the robustness of the CIs' thresholds are worse than the Li's. On the other hand, ratios below one 
mean that the CIs' thresholds work better in the experiment than the ones associated with the LI. The 
variational CIs were introduced to improve the performances of the LI. Therefore, we expect ratios below 
one, and this is the case for many of the CIs. 

The normalized time evolution of the difference between the percentages of chaotic orbits according to 
a change in the thresholds by ±1% is shown on the top panels of Fig. [TJ The weakness of the theoretical 
fixed threshold for the MEGNO becomes evident (green linepoints on the top left panel of Fig. [TJ. This 
weakness is a logical consequence of the asymptotic nature of the threshold. On the other hand, we find 
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that the RLI and the SALI have very rehable thresholds (i.e. ratios below 1), despite their empirical nature 
(red linepoints on the top left and top right panels of Fig. [T]for the RLI and the SALI, respectively), and 
the FLI/OFLI and the GALI4, as well (blue linepoints on the top left panel of Fig. [T]for the FLI/OFLI, 
and green, blue and violet linepoints on the top right panel of Fig. [T]for the GALI2, the GALI3 and the 
GALI4, respectively). The GALI4 starts with a value higher than one, but after a short initial transient, 
its rate decreases below one. 

Now we consider the percentages of chaotic orbits. The value of a final percentage of chaotic orbits 
is not as important as the approximation rate to that value. The former can be fixed by an empirical 
adjustment. The latter is a combination of the threshold's reliability and the efficiency of the indicator. 
Therefore, it is not easy to adjust in order to improve the results. Notwithstanding the relative importance 
of the CIs' final percentages of chaotic orbits, we would like to test the reliability of the thresholds given in 
Tabled) We consider the percentage of chaotic orbits given by the LI by 10^ u.t. (i.e. ~ 39.92%) the "true 
percentage" of chaotic orbits. We support this idea because the LI is the most tested CI in the literature 
and 10'* u.t. is a confident convergent time for this indicator. 

To take these considerations into account we are going to normalize the time evolution of the percentage 
of chaotic orbits according to every method with the "true percentage" given by the LI. Thus, values higher 
than 1 show percentages of chaotic orbits higher than the "true percentage" given by a confident value of 
the LI. 

The time evolution of the percentage of chaotic orbits according to the thresholds of Table [1] normalized 
by ~ 39.92% is shown on the bottom panels of Fig. [TJ The disagreement with the Li's final percentage 
of chaotic orbits (see the bottom left panel of Fig. [1]) forces an imme diate empirical adjustment of the 
MEGNO's threshold. Actually, on considering the recommendations of Maffione et ad . l201ll |. a threshold 
value close to 3 seems to be more appropriate. The SALI has the robustest threshold according to the 
top right panel of Fig. [1] (red linepoints), as we have just seen. However, the convergency of both SALI 
and GALl4's to the Li's final percentage of chaotic orbits (see the red and violet linepoints on the bottom 
right panel of Fig. [T]for the SALI and the GALI4, respectively) is slower than the convergency of the RLI 
and the FLI/OFLI (see the red and blue linepoints on the bottom left panel of Fig. [T]for the RLI and the 
FLI/OFLI, respectively). The final percentages of the RLI and the FLI/OFLI are higher than that of the 
LI, though. Nevertheless, the final percentage of chaotic orbits can be fixed with a small adjustment of 
their thresholds, as previously mentioned. 

Thus, among the above-mentioned CIs, the RLI and the FLI/OFLI show the most reliable thresholds 
in the experiment and, on comparison, the FLI/OFLI's threshold seems to work better. 

By 10^ u.t. the threshold taken for the GALI4 (i.e. t~^) reaches the computer's precision of the computer 
(10^^^) and thus, every chaotic orbit lies beyond such precision. Therefore, the last point on the right panels 
of Fig. [1] (violet linepoints) falls apart from the tendency. 



3.1.2. Qualitative study of a sample of orbits through the CIs' final values 

In this section we only use the information provided by the CIs' final values and times of saturation to test 
their speed of convergence and resolving power. 

The FLI and the OFLI, the SALI and the GALIs increase or decrease exponentially fast for chaotic 
motion (see Section [2.21 for the FLI/OFLI and Section [2.41 for the SALI/GALIs, and references therein). If 
the chaotic nature of an orbit is well characterized when the CI reaches a certain value, it is worthless to 
continue the calculation, and the integration process should be stopped. Such value is a saturation value and 
an immediate consequence of that saturation value are the times of saturation. The times of saturation are 
the time by which the CIs' final values reach the corresponding saturation va l ue. Th e tim es of saturation 



recove r the hyperbolicity levels of the chaotic component (see [Skokos et al\ . |2007I | and [Maffione et al. 
2nilal ]l. 



The parameters used for the computation of the CIs and the thresholds (see Tabled]) are the same as 
those applied in Section [3.1.11 We present the performances of the CIs on the sample of 1.25751 x 10^ i.e. 
on the energy surface defined by /i = 0.118 constraining the values of x, y and to x = 0, y G [—0.1 : 0.1] 
and py G [—0.05 : 0.05] (the same used in Section [3.1.ip with an integration time of 10^ u.t. 
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Figure 1. The normalized time evolution of the difference between the percentages of chaotic orbits according to a change 
in the thresholds of Table [1] by ±1% (a) for the RLI, the MEGNO and the FLI/OFLI, and (b) for the SALI, the GALI2, 
the GALI3 and the GALI4. The time evolution of the percentage of chaotic orbits according to the thresholds of Table [1] 
normalized by ~ 39.92% (c) for the RLI, the MEGNO and the FLI/OFLI and (d) for the SALI, the GALI2, the GALI3 and 
the GALI4. See text for further details. 



Figure [2] shows the performances of the LI and the RLI (top left and top right panels, respectively), 
the MEGNO and the D (middle left and middle right panels, respectively) and the OFLI's final values 
and times of saturation (bottom left and bottom right panels, respectively). The left panels of Fig. [3] show 
the phase space portraits of the GALI method for the above-mentioned region. The right panels of Fig. [3] 
show the corresponding times of saturation (from top to bottom panels, the GALI2, the GALI3 and the 
GALI4, respectively). 

The FLI and the SALI have not been included because they have similar performances to those of the 
OFLI and the GALI2, respectively. 

First, we deal with the identification of two domains in the chaotic component, namely the chaotic 
sea for y < —0.05 and the stochastic layers surrounding the main stability islands. We show that the CIs 
can be divided into two groups according to the different ways of describing such a component. In the first 
group we include the LI, the MEGNO and the D; in the second group, the RLI, the OFLI and the GALIs. 

Both the LI (top left panel of Fig. [2]) and the MEGNO (middle left panel of Fig. [2]) identify the two 
domains with their final values. However, the MEGNO gives a clear distinction between them. The D 
(middle right panel of Fig. [2]) does not show such a clear distinction between both chaotic domains, but it 
shows a structure in the chaotic sea which is not reflected in the portraits of the LI or the MEGNO. 

In the second group, the RLI (top right panel of Fig. [2]) shows a slight difference between the colours 
used to identify the two domains due to its high speed of convergence for chaotic motion. The OFLI 
(bottom left panel of Fig. [2]) and the GALIs (left panels of Fig. [3]) make no distinction at all, because 
of the same high speed of convergence for chaotic orbits. However, the OFLI and the GALIs count with 
the times of saturation to recover the hyperbolicity levels (bottom right panel of Fig. [2] for the OFLI and 
right panels of Fig. [3] for the GALIs). Thus, the distinction between both chaotic domains is perfectly seen. 
Moreover, the structure seen by the D in the chaotic sea is also recovered. 

The OFLI and the GALIs show the best way to describe the chaotic component thanks to the com- 
bination of a high speed of convergence and their times of saturation. Nevertheless, a similar quantity to 
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the times of saturation could be defined for the RLI. 



(a) (b) 
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Figure 2. Phase space portraits with 1.25751 x 10^ i.e. in the region defined by a; = 0, j/ £ [—0.1 : 0.1] and py £ [—0.05 : 0.05] 
in the HH potential on an energy surface h = 0.118. The phase space portraits are computed by 10^ u.t., and the color scale 
represents the final value of each indicator, (a) for the LI, (b) the RLI, (c) the MEGNO, (d) the D and the OFLI. In case of 
the OFLI we include (e) the final values and (f) the corresponding times of saturation by a final time of integration of 10^ 
u.t. in logartithmic scale. 



Second, we deal with the description of the regular component and show that the CIs can be divided 
into two groups again. The division comes from the appearance of spurious structures, which might be 
due to different reasons: numerical artifact s, short final integr ation times and a poor selection of initial 
conditions of the variational equations (see [ Barrio et al . 20091 ] for a thorough discussion). 

The LI, the RLI and the D show spurious structures in the regular component (Fig. [2]). Those structures 
are also seen in the GALIs' phase space portraits (Fig. These structures are seen as lines cutting the 
islands of stability in half. 

On the other hand, neither the MEGNO (because of its universal value for quasiperiodic motion, ~ 2, 
it does not show any structure for the regular component) nor the OFLI presented spurious structures in 
the experiment. 

Finally, we deal with the separation of both chaotic and regular components and the speed of conver- 
gence. The distinction between domains of chaotic and regular orbits is successfully drawn by every single 
method of the package. However, the OFLI and the GALIs seem to make the distinction faster. 
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Figure 3. Left panels: phase space portraits of the GALl method with 1.25751 x 10^ i.e. in the region defined by x = 0, 
y G [—0.1 : 0.1] and py € [—0.05 : 0.05] in the HH potential on an energy surface h = 0.118. (a), (c), (e) for the phase 
space portraits by 10^ u.t. of the GALI2, the GALI3 and the GALI4, respectively, (b), (d), (f) for the corresponding times of 
saturation by a final time of integration of 10*^ u.t. in logarithmic scale. The color scale has the same meaning as in Fig. [2] 

In the case of the GALI^ , as we increase the index k of the method (from top to bottom panels of Fig. 
[3]), we find more detailed descriptions of the regular and chaotic components by means of the GALIs' final 
values and the corresponding times of saturation, respectively. Furthermore, the GALI4 has the highest 
speed of convergence in the chaotic orbits. 

Therefore, a combination of the FLI/OFLI and the GALI4 seems to carry enough information to 
successfully describe the strongly divided phase space of the studied region in the HH potential by means 
of their final values and times of saturation. 



3.2. Study of some representative orbits through the CIs' time evolution 

In [ Maffione et"ZI . l2nilal ]. Section 3, the authors used the CIs' time evolution to study the spe ed of con- 



vergen ce and the sensitivity to distinguish different kinds of motions. The results found by [Maffione et at 



2011al ] for mappings remain the same for the HH potential. Therefore, we deal with some CIs' features not 



discussed in their work. We start studying the ability of the CIs to distinguish periodic orbits and their 
CIs' dependency on the i.d.v. Then, we deal with the sensitivity of the GALIs to identify instabilities. 
We study a small sample of six orbits in the HH potential: 
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Orbit (a): the 5-periodic orbit with i.e.: ~ (0,0.35207,0,0.14979), and two close 

orbits on the same energy surface, h = 0.125. The close orbits are defined by Ay = 0,00793 and Ay = 

0. 02793, Ay being the difi^erence in the y-con iponent of the i.e. of orbit (a) and the y-component of the 

1. e. of the close orbits. They were taken from [Manos et al . 201 1[ |. 

Orbit (b): a quasiperiodic orbit with i.e.: {x^^''\ y^^^\ px^^\ py^''^ ) ~ (0,0,0.5,0), which is associated with 
the 1-periodic orbi t in the HH potential on the energy surface h = 0.125. The orbit was taken from 



Skokos et a/ll2007l ]. 



Orbit (c): a regular orbit close to the separatrix with i.e.: 
Orbit (d): a chaotic orbit inside the chaotic sea with i.e.: (x^^'^^ 



(0,0.5085,0.25512,0). 
(0,0.6,0.14142,0). 



The last two orbits belong to an energy surface of /i = 0.118 and were taken from [Cincotta k, SimdI . 

1999l |. 

The final integration time is 2.4 x 10^ u.t. (unless stated otherwise) and ensures the eonvergeney of 
the LI for both energy surfaces to properly determine the regular or chaotic nature of the orbits. 



3.2.1. Identification of periodic orbits 

We start with the identification of periodic orbits. Thus, we evaluate the CIs in the orbit (a) and the 
quasiperiodic orbits close to (a). The behaviors of the LI, the SSNs, the MEGNO, the FLI, the GALI3 
and the GALI4 are nearly the same for the three orbits. 

The OFLI is precisely defined to reveal periodicity ( jPouchard et al. . 2002l |). The top left panel of 
Fig. m shows the performances of the OFLI for the orbit (a) and for the two quasiperiodic orbits close to 
(a). Although, the OFLI shows an oscillatory regime around a constant value for the periodic orbit, an 
increment following a linear rate is seen for quasiperiodic motion, as in the case of the FLI. 

The top right panel of Fig. H] shows the performances of the GALI2 (as the SALI has similar behavior, 
it is not included) on the periodic and quasiperiodic orbits of the sample. The GALI2 is the only GALI 
capable of clearly distinguishing periodic motion in the experiment. 

Let us analyze the behavior of the GALIs (and the SALI) for the special case of a 2-d.o.f. Hamiltonian 
system as the HH potential. In a 2~d.o.f. system, the quasiperiodic orbits lie on 2-dimensional tori. All 
deviation vectors tend to fall on the 2-dimensional manifold tangent to the torus where the motion takes 
place. If we start the computation with any two linearly independent deviation vectors (in order to compute 
the GALI2 or the SALI), they will remain linearly independent on the 2-dimensional tangent space of the 
torus. Then, the GALI2 (or the SALI) is constant and different from zero. On the ot her hand, the GALI^ 
(with /c = 3, 4) tends to zero, since some i.d.v. become linearly dependent on the rest i; [Manos et oll . l201l[ |). 

In the case of the periodic orbits, the motion takes place on 1-dimensional tori (an invariant curve), 
the tangent space being also 1-dimensional. Thus, the behavior of the GALI2 (or the SALI) is oc t~^. The 
following asymptotic formulae summarizes the general behaviors in the HH potential of the GALI^ for m 
initially tangent i.d.v. and regular orbits lying on M-dimensional tori (with M = 1,2): 



' constant i/ A; = 2, M = 2, m = 0, 1, 2 



GALIk{t) oc < 



ifk = 
ifk = 

ifk 
ifk 



t ^ ort 



2, M = 

3, M = 

= 3,M 
= 4,M 

4,M 



1, m = 

2, m = 

= 2, m 
= 2, m 



0,1 
1,2 

= 
= 2 



i/fc = 4,M = 2,m = 0,1 

For the general formulae see j Christodoulidi Bountid . l200fil | or Section [231 of this paper. 
On the top right panel of Fig. [U we see that the GALI2 shows the same initial behavior for the three 
orbits (i.e. following a power law of which is also included in the figure). After an initial transient of 
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similar behaviors, the GALI2 oscillates around a constant value for the quasiperiodic orbits. 

In a 2-d.o.f. Hamiltonian system, the D is not able to distinguish between periodic and chaotic mo- 
tion (a decay is observed for both moti ons). N e verth eless, the ambiguity disappears when dealing with 
Hamiltonian systems of more d.o.f. (see [Skokosl . I2OOII ] for further information). The bottom left panel of 
Fig. U shows the behavior of the D for the two regular orbits close to the orbit (a) and for an extended 
integration time of 10^ u.t. The closer the orbit to the periodic orbit, the longer it takes the D to converge 
to a constant value, characterizing the regularity of the orbit. 

On the bottom right panel of Fig. Owe show the performances of the RLI. The RLI converges to lower 
values for the orbit (a) than for the quasiperiodic orbits. The closer the quasiperiodic orbit to the orbit 
(a), the lower the final value of the RLI. Nevertheless, neither is there a significant change in the behavior 
of the CI nor a reference value is defined to separate the periodic orbits from the associated quasiperiodic 
orbits. 
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Figure 4. Performances of the CIs for the orbit (a) and for the quasiperiodic orbits close to (a), (i) for the OFLI and (ii) 
the GALI2 (the predicted behavior for the periodic orbit is also included in the left chart). Performaces of the CIs for the 
quasiperiodic orbits close to the orbit (a), [Hi) for D and [iv) the RLI. 



According to the previous results, the OFLI and the GALI2 (and the SALI) are the only CIs capa- 
ble of clearly identifying the periodic orbit. The identification of the periodic orbits may fail due to an 
inappropriate c hoice of the i.d.v. Thoug h the dependency of some fast CIs on the i.d.v. has already been 
pointed out by (Froeschle &: Legal . 200Cll | and (Bar rid . 2005], it is worth reviewing this dependency. 



3.2.2. Dependency of various CIs on the i.d.v. 

In Section [3.11 we used the canonical basis of i.d.v. for the whole sample of orbits along the experiments. 
However, herein we analyze only the orbit (b) and we use three different bases of i.d.v. to test the de- 
pendency of the CIs on the i.d.v. These bases are characterized by the parameter m (the number of i.d.v. 
initially tangent to the torus). The three bases are selected to have m = 0, 1 and 2 (2 is the upper limit 
in a 2-d.o.f. Hamiltonian system like the HH potential). From [Skokos et al. . 200?! ] we know that the unit 



vectors (1,0,0,0) and (0,0,0,1) are initially tangent to the torus where the orbit (b) moves. Then, using a 
Gram-Schmidt process we build three bases of i.d.v. with m = 0, 1, 2. For example, the basis with m = 2 
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has four orthonormal vectors, two of them being the above-mentfoned i.d.v. 

Although the LI and the RLI show some dependency on the parameter m, the differences are mean- 
ingless. The MEGNO, the FLI and the OFLI do not show any dependency at all. 

The top left panel of Fig. [5] shows the behavior of the D for the orbit (b) and for the three bases 
of i.d.v. We identify the regular nature of the orbit (b) using the bases with m = 1,2. The D reaches a 
constant value with both bases, with different convergent times, though. On the other hand, with m = 0, 
the D decreases like for chaotic motion. Despite extending the integration time to 10^ u.t., the result has 
not changed. 

To study th e beha viors of the GALI^ with A; = 2, 3, 4, we reproduced the experiment shown in Fig. 4 
of [Sk okos et ai . 2007]. They analyzed the orbit (b) in the HH potential with m = 0, 1,2. Although the 
results agree with m = 1, 2, we see a discrepancy in the behavior of the GALI2 with m = (see the dark 
gr ay line on the top ri ght panel of Fig. [5] and their top left panel of Fig. 4). On the top left panel of Fig. 4 



m 



Skokos et al 



2007l |. the authors showed an almost constant value for the SALI and the GALI2 (we do 
not include the SALI in our figure because the results are similar to the GALI2). Although we select zero 
initially tangent i.d.v. as they did, our two i.d.v. used for the computation of the GALI2 align with each 
other. On the bottom left panel of Fig. [5] we show the opposite behaviors of the quantity cos(a) (a being 
the angle between the direction of the i.d.v. and the direction of the flow) for both i.d.v. Therefore, they 
are linearly dependent and we see a power law for the GALI2 (top right panel of Fig. [5]), like for a 
periodic orbit. This linear dependency on both i.d.v. also explains the behavior of the D for the orbit (b) 
and with m = 0. 

On the bottom right panel of Fig. [5l we show the behaviors of the quantity cos(a) for the i.d.v. (1,0,0,0) 
and (0,0,0,1). They correspond to the orthonormal basis of i.d.v. with m = 2. It is clear that the first i.d.v. 
aligns parallel to the flow while the second i.d.v. oscillates around one of the orthogonal directions to the 
flow. The last case wit h m = 2 is the general case, an d th e behaviors of the GALIs are well predicted by 
the formulae given by Christodoulidi &: Bountisl . 2006l | or Skokos et al. . 2007]. Therefore, there is no need 
to check the choice of the i.d.v. to analyze large samples of orbits. It is convenient to track the behavior of 
the i.d.v. for a short period of time to analyze small sample of orbits, though. 




Figure 5. (?) Performaces of the D for the regular orbit (b) and m = 0, 1, 2. (ii) GALIs' performances for the orbit (b) with 
m = 0. The power laws associated with each GALI are also included (note the logarithmic scale). Behavior of the quantity 
cos(a) for both the i.d.v. considered in the computation of the GALI2 (m) with m — and (iv) with m — 2. 
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3.2.3. Sensitivity of the GALIs on dynamic instabilities 

We are also interested in studying the CIs' sensitivity to instabilities. In Maffione et al . 2011al ] the authors 
dealt with the sensitivity to stability and chaoticity levels of many of the CIs of the package, i.e. the LI, 
the RLI, the MEGNO, the SSNs and the D, the FLI and the SALI. Although their study was done on 
a 4D mapping, the results do not change in the HH potential. Furthemore, the FLI and the OFLI show 
no difference in their performances. Therefore, here we deal with the sensitivity to the instabilities of the 
GALIs alone and, for this purpose, we take the orbits (c) and (d). 

On Fig. [6] we show the behaviors of the GALI2, the GALI3 and the GALI4 for the orbits (c) (left panel) 
and (d) (right panel). The GALIs for the orbit (c) show an unstable behaviour for a regular orbit because 
of the large amplitude of the oscillations around the predicted behavior for regular motion (i.e. GALI2 oc 
constant; GALI3 oc t~^; GALI4 oc t~^, with m = 2). On the other hand, the GALIs for the orbit (d) 
show a chaotic orbit because their behavior follows exponential laws (i.e. GALI2 oc e~^^^; GALI3 oc e~^-^^*; 
GALI4 oc e-^^i*, with xi ~ 0.4118490 x 10"\ xi being the LI value of the orbit). The GALI4 has the 
highest speed of convergence for the chaotic orbits among all the indic ators tested so far. However, it is 
not economical in terms of computational time, see Skokos et al\ . l2007l | and Section 13.31 of this paper for 
further details. 




Figure 6. (i) Performances of the GALIs for the orbit (c). (u) Performances of the GALIs for the orbit (d). The predicted 
behavior is also included, in logarithmic scale, in both charts. 



Finally, we show that the OFLI and the GALI2 (and the SALI) are the only CIs that clearly identify the 
periodic orbit (a) due to an evident change in their behaviours (Section l3.2.ip . However, some dependency 
on the i.d.v. it is shown in the experiment, which may lead to an incorrect identification of periodic orbits 
(Section I3.2.2p . Thus, it is convenient to track their behaviors for a short period of time in order to select 
the bes t i.d.y . The dependency of first order variational indicators on the i.d.v. is further discussed in 
|Barriol . [2005l |. 

The GALI2, like the other GALIs, shows a dependency on the param eter m which gives information 
of the dimensionality of the tori where the regular motion takes place (see Manos et al. . 201 1[ |). Moreover, 
in the case of the chaotic orbits, the GALI4 has the highest speed of convergence among all the indicators 
tested so far (Section I3.2.3p . 

The MEGNO showed the same good performan ces for single studies of orbits in the HH potential (not 
shown) as for mappings (see [Maffione et adl2011al ]). in particular because it is simple to characterize the 



stability levels, see Section [2.11 Thus, for studies of this kind, a combination of the OFLI and the MEGNO 
seems to be the best choice if we can use the time evolution of the CIs. However, the GALI2 (or the SALI) 
is also appropriate for the task. 



3.3. The computing times 

Before going to the second part of the work, we finish the comparative study of variational CIs studying 
their computing times. 
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On Fig. [71 we show the CPU-times for ah the CIs currently under study. On the left, we compute the 
CIs for a sample of 10^ i.e. located in the interval defined by y € [0.55 : 0.6], x = py = Q and h = 0.118 
(hereafter region A). Region A has a divided phase space, and a high percentage of orbits that reach the 
saturation values (i.e. 10^° for the FLI/OFLI and lO"^*^, the precision of the computer, for the SALI or the 
GALIs) within the interval of integration (1.2 x 10^ u.t.). On the right, we computed the CIs for a sample 
of 3 X 10^ i.e. located inside the stability island associated with a periodic orbit of period one. We call this 
region "B"; it is defined by y G [0.3 : 0.506], x = py = and h = 0.118. In this region there are no orbits 
that reach any of the saturation values by 1.2 x 10"^ u.t. 

From Fig. [7] it is clear that the CPU-time increases with a growing percentage of regular orbits for 
each CI (compare both panels). Also, the CPU-time increases with the complexity of the CPs algorithm. 
Therefore, the LI is the least time consuming method while the D is the most CPU-time demanding CI. 
The order would be: the LI, the MEGNO, the SALI and the RLI, the FLI and the OFLI, the GALI2, the 
GALI3, the GALI4 and the D. The MEGNO uses two more equations than the LI. The SALI and the RLI 
have very similar computing times because the SALI uses the computation of two deviation vectors while 
the RLI, the computation of two orbits. The FLI and the OFLI use the computation of four deviation 
vectors (the HH potential has a 4D phase space). The GALIfc uses k deviation vectors, but their numerical 
implementation is time consuining (d espite the fact that we used the "Singular Value Decomposition" or 
SVD routine, see Skokos et al. . 20081 ]). Finally, the D uses a statistical approach which is not numerically 
efficient. 




Figure 7. CPU-times for all the CIs under study (a) for region A and (b) for region B. 



Using a saturation value for some methods has several advantages which were studied in this work 
(Section 13. 1.2p . like recovering the chaoticity levels with the times of saturation. Furthermore, there is also 
a numerical advantage with such kind of methods. If the orbit r eaches the saturatio n valu e of the CI, the 
compu tation can be stopped, thus reducing the CPU-time (see [ Skokos et al. . 20071 ] and Maffione et al. 



2011al ]). On Fig. El we tried to measure such relationship. Therefore, we computed the CPU-time gained 



due to the percentages of orbits that reach the corresponding saturation value. This relation gives an idea 
of how much efficient an indicator is. For example, with a 35% of orbits that reach the saturation value 
10"^*^, the SALI gains less than 5% of the CPU-time, the GALI2 and the GALI3 gain around 5% and 12% 



0, 



respectively and the FLI and the OFLI gaining between 20% and 25%. Although the GALI4 shows higher 
percentages of orbits that reach the associated saturation value (10^^®), the CPU-time gained is not as 
high as that for the FLI or the OFLI. 

Finally, this experiment shows that the FLI/OFLI is the most efficient method because it gains most 
computing time by means of a succesful implementation of the saturation values. 
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Figure 8. Percentage of orbits that reach the saturation values against the reduction in the CPU-times for all the CIs that 
have a saturation value. 



4. The software package using Taylor method 

Here we summarize the main features of the software p ackage that integra tes ODEs using the Taylor 
method (tayloiH hereafter) that is described in detail in [jorba fc Zou . 2005l |. 

As mentioned in Section [H one of the main distinguishable features of this software is the use of 
the so-called automatic differentiation^ which computes the derivatives of the functions involved up to 
an arbitrarily high order. This particular implementation of automatic differentiation generates its code 
from a file with the differential equations system, and makes the differentiation in an automatic fashion, 
without losing accuracy on the d erivatives. There a r e many general pur poses computer prograrn s which 
build its code from a program f [Beda et all Il959l |. jCibbond . Il960l | & IChan g fc Corlisd . Il994l |l. while 
others are built up by the user for a specific problem (see for instance, iBroucke . Il971| for the A^— body 
pr oblem). Further alt e rnativ es to apply the Taylor method can be found in jSavageau &: Voit . Il987| | and 
in [irvin e &: Savageau . 199[lf |. In fact, the key point of the Taylor method is the precise computation of the 
derivatives of x{tm), x^^^tm), which the automatic differentiation plays a crucial role. 

The way in which t ay lor operates can be summarized as follows: the package reads the differential 
equations system from a plain text file, and uses it to build the integrator (for a full syntax and parameters 
list see "Taylor user's manual", which comes with the package taylor). taylor subroutines are written in 
ANSI C code. It includes the subroutines for the automatic differentiation, the time stepper, as well as step- 
size and order control. The subroutine which computes the automatic differentiation reads each equation 
from the file, decomposes the expression into binary or unitary operations and gen erates recurre r ices i n 
order to compute the derivatives analytically (for an example of some recurrences see |jorba fc Zoul . l2005l |. 
Section 2, Proposition 2.1). Thus, the derivatives of the function / up to an arbitrarily high order can 
be computed without losing acc uracy. For a more detailed description of the software implementation see 
Section 4 of j.Torba fc Zoi] . l2nn,^ . 



For the calculation of the step size h and the order of the method p, taylor uses variable-stepsize 
and variable-order algorithms. The package computes, for each time-step, the order pm that guarantees a 
truncation error of the order of a given threshold £m (for simplicity, in the numerical computation sections 
we call it e). With this value of pm the step size hm is computed. These two values are chosen in such a 
way that not only the required precision e is achieved but the number of operations needed to advance one 
time step is minimized. 

An interesting feature of Taylor method, unlike other integrators, is that to increase accuracy, it 



*The package taylor is released under the GNU Public License, and it can be retrieved from 
http : // www . ma . utexas . edu/ ~mzou/ taylor (US) 
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re quires much le s s wor k to increase the order p than to reduce the step size h. Its explanation can be seen 
m jjorba fc Zoul . l2005l ]. on Section 3.3. 



Another useful feature of taylor is the possibility to implement easily the use of extended arithmetic, 
using high precision libraries like GMP^. This is possible due to the implementation of abstract operations 
with a floating point type called MY_FLOAT. Then, in the header file, these operations are replaced by the 
ones corresponding to the mathematical library we are using (the default is double). 

5. Application to a model of triaxial elliptical galaxy 

We considere d a t riaxial potential, corresponding to an auto-consistent model of elliptical galaxy 



Muzzio et al\ . l2005l | , given by 



= -f,{p)-f,{p){x' - /) -Up){z' - y% (7) 

where p is the softened radius 

/ = r^+6f for/o, (8) 

and 

p^ = r^ + ^ef for/,,/,, (9) 

with the softening parameter fixed at ei = 0.01. 

The functions fs were computed from an A^— body simulation and fit with equations of the form 

C 

fs = -n: s = 0,x,z. (10) 

(p " + qg^y'l''^ 
The adopted values for C^, kg, Qs, Is are given in Table 2 

Table 2. Values of the coefficients of the functions 
fs{p) given by Eg. [Tol (taken from Muzzio et ai, [2005]). 







Cs 


ks 


qs 


Is 


s = 





0.92012657 


1.15 


0.1340 


1.03766579 


s = 


X 


0.08526504 


0.97 


0.1283 


4.61571581 


s = 


z 


-0.05871011 


1.05 


0.1239 


4.42030943 



We adopted value E = —0.5 for the system's energy (which corresponds to an x-axis central orbit 
with period ~ 10 u.t.). We considered central orbits (i.e. xq = uq = zq = 0) and, by defining AV = 
E — V{xo,yo, zq), we took initial values for px, Py and pz which satisfy px'^ + Py^ + Pz^ = 2Ay using a 
grid made of points which swept the values of px and py with a step of Apx = Apy = 0.2, yielding a 
total of 2080 orbits. The integrations were carried out for t = 5 x 10^ and t = 5 x 10^ u.t., but only the 
results for t final = 5 x 10^ are shown here since the same behavior is observed in both cases. We chose as 
Unitiai = 0.1 in Order to avoid a division by zero in the computation of the MEGNO (note that there is a 
t in the denominators of eqs. [4]and[5]), and tinuiai = in Subsections 15.21 and 15.31 The required tolerances 
were e = 10~^^ (the same for both the absolute and relative tolerances). 

For all the computations done in this second pa rt of the work we us ed variable-stepsize. For taylor, 
we used the stepsize control method mentioned on [Jorba fc Zot1|200,5I ]. Section 3.2.1, and for D0PRI8 
and BS their own stepsize control method, included on both programs. All the results presented here have 
been computed in a 2xDual XEON 5120, Dual Core 1.86 GHz. 



^For further informatin see http://gmplib.org/ 
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5.1. The MEGNO 



The skillful computation of the MEGNO definitely requires a n efficient integrator . Therefore, it would be in- 
teresting to confirm the superiority of taylor claimed out in j.Torba fc Zoul . l2nn,^ . To this end, the package 
was compared with o ther t wo well-known integrators: a Runge-Kut ta 7/8 method, the so-called D0PRI8 
JPririce fc Dormand . 1981], and the Bulirsch-Stoer method (BS) [Numerical Recipes in FORTR AN TtI . 
' l997l |. The results of the comparison regarding both speed and precision of the integrations performed 



by means of taylor, D0PRI8 and BS are presented in Figure [U where we have plotted the respective 
CPU times and maximum energy errors for one out of ten i.e.. 




1000 

No. of initial condition 



1000 1500 
No. of initial condition 



Figure 9. CPU times (Left) and Maximum energy error (Right) for several i.e. in the calculation of the MEGNO for the 
triajdal potential, using taylor, D0PRI8 and BS, with tjj_nal = 5 x 10^ u.t. 



The left side of Figure [9] reveals that taylor is considerably slower than D0PRI8 and BS. Its computation 
time is ~ 2 or 2.2 times the time insumed by the other two integrators (i.e. it's slower in a factor of ~ 2, 
2.2). When we look at the maximum energy error displayed on the right side of Figure [9l we observe that 
taylor is slightly more precise than D0PRI8 (in a factor ~ 2), while the errors arising from the integration 
by BS are about 2 orders of magnitude larger (in a factor ~ 70). 

Therefore we can assert that although taylor is nearly as precise as D0PRI8, it clearly demands much 
more computation time than both D0PRI8 and BS. Thus, at least for this particular system of ODEs, the 
taylor package turns out to be a somewhat inefficient integrator. 

This is a surprising fact, since taylor showed up as a very competitive integrator in [jorb a fc Zou , 

20051 ]. However it seems that there are some particular features in these ODEs that could be conflictive for 
the taylor package. In fact, the differential equations for the MEGNO 



{7/ _ j. h{xo;t) 
-^<57(a,o;t) (11) 

encompass the temporal variable t in the expressions for Z' and Y' in the numerator and denominator, 
repectively. As a consequence, we infer that, at the beginning of the integration when t is very small and 
near its end, when we consider a large period, we obtain large values of Z' or Y' , respectively, taylor 
might compensate this effect by taking very small step sizes, and these two equations might be responsible 
for the inefficient performance of the package. 

In order to test this conjecture we proceed to the computation of another dynamical indicator, namely, 
the LI. This suits perfectly our needs, since only the equations of motion and their first variational ones 
are to be integrated, and the two equations (fTTIl that calculate the MEGNO are no longer involved. 
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5.2. The LI 

The task was carried out using the very same i.e. and parameters for the galactic model as those of the 
previous section and the same tolerance was imposed. The CPU times demanded by the computation of 
the LI values and the resulting maximum energy errors are displayed in Fig. [TUl 




Figure 10. CPU times (Left) and maximum error on the energy (Right) for several i.e. in the calculation of the LI for the 
triaxial potential, using taylor, D0PRI8 and BS, with t^i„(j; = 5 x IC^ u.t. 

Regarding to the CPU time, taylor still remains slower than both D0PRI8 (a ~ 1.4 factor) and BS 
(~ 1.65), but it is evident that the differences have decreased significantly, i.e. from a ~ 2 to ~ 1.4 factor. 
However, as long as the energy error is concerned, no significant dissimilarity from the results shown in 
Figure [9] is observed. 

These results support our assumption that the equations for the MEGNO decrease taylor's perfor- 
mance. In fact, the reduction of the system of ODEs by avoiding those required for the MEGNO con- 
siderably dimished the difference between the computational times invested by the compar ed integrators . 



None theless, taylor kept on being the least efficient, in discrepancy with the results given in [Jorba &: Zou . 
imi, Table 4. 

The question arises whether for integrating only the system o f equations of motio n, taylor succeds in 



showing its benefit regarding the speed of computation stated by [Jorba &: Zoul . l2005l | . Table 4. This issue 
is addresed in the forthcoming subsection. 

5.3. Integration of the equations of motion 

On the left of figure [TT] it is displayed the CPU time insumed by taylor, D0PRI8 and BS for several i.e. 
when we integrate the equations of motion for the galactic model. The plot evinces that for this particular 
problem, taylor results slightly faster than both D0PRI8 (a factor of ~ 1.06) and BS (~ 1.08). 

Regarding to the maximum error in the energy, a similar behavior is revealed on the right side of Fig. 
[TTt being taylor slightly more accurate than D0PRI8 and about 2 orders of magnitude more precise than 
BS. 

Let us remark that when solving the three different sets of ODEs for this particular system (the 
MEGNO, the LI and the equations of motion), taylor showed up to be as precise as D0PRI8 and much 
more precise than BS; but as far as the required CPU times are concerned, only when integrating the 
equations of motion taylor turned out to be the best, in a smaller degree, though. 

The results achieved so far suggest that taylor is not well suited for dealing with systems involving 
variational equations or ODEs like the ones defining the MEGNO, but it offers some advantage when only 
the equations of motion are integrated. 

In fact, the taylor package consumes a quite significant amount of CPU time to compute the automatic 
differentiation, build the power series and evaluate them at the requested points. Thus, it seems fair to 
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No. of initial condition No. of initial condition 

Figure 11. CPU times (Left) and maximum error on the energy (Right) for several i.e. when computing the equations of 
motion for the triaxial potential, using taylor, D0PRI8 and BS, at tfi^al = 5 x 10** u.t. 

wonder whether using this dynamical system for the comparison is a fortunate selection, due to the rather 
cumbersome expressions involved. 

6. Application to a simpler problem: the perturbed 3D quartic oscillator 

In this section we present the results of the reproduction of the same numerical experiments for a somewhat 
simple dynamical problem, namely, a 3D quartic oscillator subject to a cubic perturbation of the form 
ax^{y + z), where a = 5 x 10~^ is the perturbative parameter. Therefore, the potential can be recasted as 

V{x, y, z) = + / + z^) + ax\y + z). (12) 

For these experiments we adopted the energy value E ~ 0.485, for which the central x-axis orbit has 
a period of 27r u.t., and we constructed a lattice of i.e. by taking y = Px = Pz = and values of x and z in 
the range < x < 1.5 and —1.5 < z < 1.5, with step-sizes Ax = Az = 0.05, yielding py^ > (i.e. py G M). 
The integrations for the resulting 1056 i.e. extended over lO'^ u.t. and 10^ u.t., though only the results for 
t = 10^ are included herein, for the same reason as in the case of the t final for the triaxial potential. As 
in Section [5.H the initial time tinuiai = 0.1 was adopted to avoid a division by zero in the calculation of 
the MEGNO, but the value tinUiai = was fixed in Sections 16.21 and 16.31 The tolerance was again set to 
£ = 10-1^ 

6.1. The MEGNO 

Figure [12] (left) illustrates the computational cost for evaluating the MEGNO in the current application. 
It can clearly be seen that taylor required nearly 3 times the CPU time demanded by D0PRI8 and 2.3 
times than that taken by BS for the same task. In other words, taylor turned out to be from about 2.3 
to 3 times slower than the other two integrators, a factor of the same or even larger order than the one 
observed for the triaxial potential. 

The results in Figure [T2] (right) show that, in this case, taylor is more accurate than D0PRI8 in a 
~ 15 factor, that is a difference of approximately 1 order of magnitude in the conservation of energy. On 
the other hand, the same difference as in the case of the triaxial potential (see Figure [H right) is attained 
in the current application when taylor is compared with BS performance, i.e. 2 orders of magnitude in 
favor of the former (more specifically a factor ~ 70). 

Although taylor allows better energy conservation, the time required to compute the MEGNO is too 
long to be a preferable alternative to the other two integrators, even when we are dealing with such a plain 
system as a perturbed quartic oscillator. 
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Figure 12. CPU times (left) and maximum energy error (right) for several i.e. in the computation of the MEGNO for the 
quartic oscillator, using taylor, D0PRI8 and BS, for tfi^al = 10'^ u.t. 

6.2. The LI 

The LI values were obtained for the same set of i.e. as those used in the previous section, imposing the 
same tolerance. The required computational effort and the achieved accuracy in the energy preservation 
are illustrated in Fig. [13l 
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Figure 13. CPU times (left) and maximum energy error (right) for several i.e. in the computation of the LI for the quartic 
oscillator, using taylor, DOPRI8 and BS, for tf^nal = 10^ u.t. 



From Figure [13] (left) it can be deduced that taylor was about ~ 1.2 times slower than D0PRI8 for 
the LI computation, which is a remarkably smaller difference than that observed for the CPU times for 
obtaining the MEGNO. 

Regarding the resultant maximum energy error, Fig. [T3] (right) reveals a quite different behavior from 
that observed in the case of the MEGNO evaluation. Indeed, though there is still a difference of about 2 
orders of magnitude between taylor and BS (a ~ 100 factor in favor of taylor), the difference between 
taylor and D0PRI8 grew considerably (to a ~ 45) compared to that shown on the right side of Figure 

m 

Therefore, this result differs from the one obtained for the triaxial potential, where the differences 
in energy preservation turned out to be almost the same in the three problems analysed, namely, the 
computation of the MEGNO and the LI values and the integration of the equations of motion. 
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6.3. The equations of motion 

This section is aimed at comparing both speed and precision when only the equations of motion are 
computed for the quartic oscihator by testing again the three integrators under study. 




Figure 14. CPU times (left) and maximum energy error (right) for several i.e. in the calculation of the equations of motion 
for the quartic oscillator, using taylor, D0PRI8 and BS, for t/i„a/ ~ 10^ u.t. 



As the left side of Figure [T^ shows, taylor become the fastest one. In fact, a similar behaviour with 
respect to D0PRI8 is observed as the one achieved for the galactic model (i.e. taylor is just ~ 1.08 times 
faster than D0PRI8), but the difference is greater when comparing taylor with BS in a ~ 1.3 factor. 

Concerning the maximum energy error presented on the right side of Figure [TH the differences between 
taylor and BS slightly decreased, falling below the 2 orders of magnitude, i.e. nearly a 40 factor. Moreover, 
this is the same difference between taylor and D0PRI8 results, the error curve for the latter lying behind 
the one corresponding to the BS method. 

Thus, it should be highlighted that, when we solved the three different sets of ODEs (those for the 
MEGNO, the LI and the equations of motion) for the perturbed quartic oscillator, taylor was more 
accurate than D0PRI8 and BS. In terms of CPU times, taylor was comparable to D0PRI8 and faster 
than BS when integrating solely the equations of motion, and required the largest computational time 
when determining both the MEGNO and the LI values. 

Again, the results obtained for the perturbed quartic oscillator indicate that taylor is not well suited 
for systems involving variational equations or ODEs like the ones defining the MEGNO. However, it does 
offer some advantages when only the equations of motion are integrated. The taylor performance slowdown 
is still not fully understood and it is subject to further study (extracted from a private conversation with 
Angel Jorba). 

As mentioned in Section [H the study of diffusive processes in phase space requires a fast and, above all, 
very reliable method for integrating the equations of motion of such a system. The experiments performed 
herein point out that taylor is more efficient than the other two integrators considered, at least as far 
as the integration of the equations of motion is concerned. Thus, we are bound to advance that it is the 
best candidate for studies on the accurate determination of coordinates in phase space. This issue will be 
discussed in Section [71 



6.3.1. Tolerance for a fixed AE 

In the previous section, we obtained the precision on the energy for a fixed value of e. The aim of this 
section is to obtain a value of e for which we calculate a precision of AE in the preservation of the energy. 

In order to do this we used the following steps: first, we fix the precision on the energy AEq, then we 
set an initial value for eg = AEq and, by decreasing e by one unit (i.e. 3 x 10~^^, 2 x 10~^^, 1 x 10~^^, 
9 X W-^"^, etc.), we continue iterating until we obtain a value of e for which AE < AEq. Table 3 shows 
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the results obtained. 



Table 3. Tolerance required to achieve the desired 
AS for t final = 10'' (top), tfi,^al = 10^ (bottom). 



f = 10^ 




Tolerance (e) 






AS 


taylor 


DOPRI8 


BS 




10-11 


3 X IQ-i^ 


4 X lO-i* 


4 X 10" 


la 


10-12 


3 X IQ-l^ 


4 X IQ-i^ 


2 X 10" 


14 


7 X 10-1^ 


3 X IQ-l^ 


3 X IQ-i^ 


2 X 10- 


14 


t = lO'' 




Tolerance (e) 






AS 


taylor 


DOPRI8 


BS 




10-9 


1 X IQ-ll 


3 X IQ-i^ 


5 X 10- 


12 


10-1" 


3 X 10-13 


4 X IQ-i"' 


4 X 10- 


13 


10-11 


3 X 10"13 


4 X IQ-i^ 


2 X to- 


14 



As indicated in Table 3, considerable less precision in e is required by taylor to achieve the same error 
on the energy preservation. Note that, in the case of t final = 10^, the three values of e given by taylor are 
exactly the same. T his is due to how taylor selects the order p and the stepsize h (for further details, see 
|jorba k. Zoul . 2005], Section 3.2). A brief explanation is that the order p depends on e and the stepsize h 
basically depends on p. Therefore, for all the values of e which result in the same value of p, we will obtain 
the same /i, in an identical integration. This step-like behaviour can be seen in Figure [TSl 




le-10 1e-12 1e-14 1e-16 
Tolerance (eps) 



Figure 15. Values of e for different values of AS. 

Here we can see that, for e ?a 2 x lO-i^ in D0PRI8 and e ^ IQ-i^ for BS, as both integrators saturate, 
due to their intrinsic precision, it is useless to compare them beyond those values. 

The step-like nature of taylor's curve means that, if we take a given value of /S.E, we can increase the 
tolerance up to the largest value in the same step. For instance, seeing figure [T5t if we ask iS.E ~ 3 x lO"!^, 
we can increase e up to ~ lO-i^. Or, from another perspective, if we change from e = lO-i^ to 3 x lO-i^, 
we win approximately 3 orders of magnitude, while D0PRI8 and BS will improve less than 1 order of 
magnitude. 



7. On the accuracy of the computation of coordinates in phase space 

In this section we will try to elucidate which of the three integrators under study permits the most precise 
determination of the coordinates in phase space for a given initial value problem when investigating diffusive 
processes. 
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To this end, the integration using taylor with extended precision arithmetic (i.e. GMP hbrary with 
a mantissa of 256 bits, see http://gmplib.org/ for details) could be considered as the "exact" solution 
against which to contrast the result of integrati ng the equations of motion by means of the three methods. 
This choice is based on the results obtained in j.Torba fc Zoul . l2nn.^ . Section 5.1.4. 

Of the 1056 i.e. used in Section[6l in order to perform the accuracy test, the perturbed quartic oscillator 
was integrated with those which the MEGNO classified as regular (i.e. 1039 orbits). The reason for that 
selection is that it is useless to compute the error on the coordinates of chaotic orbits, because of their 
sensibility to the i.e., resulting in too different final coordinates for each integrator. The results were 
then compared to the extended precision approximation. The integration was made for several tolerances, 
ranging from 10"^*^ to 10~^^, and stopped at 10^ u.t. 

In fact, the magnitude of the error vector in configuration space provides an idea of the precision of 
the integrator applied. Nonetheless, the magnitude of its component along the direction normal to the 
trajectory would be more relevant in the present application is, since the error in the tangential direction 
could be assimilated as a mere shift of the integration time, having no impact on a diffusion analysis. If 
we define 9 as the angle between the position error vector and the tangent line of the actual orbit, given 
by the velocity vector of the "exact" trajectory, that is. 



Ar • V = Arvcos(6'), (13) 

the tangential and transverse components of the error can be recasted as Arcos(0) and Arsin(0), respec- 
tively. 

Table 4 clearly shows the absolute value of the position error and its component in the transverse 
direction as a percentage of the whole for several tolerances. To get these results we considered all the i.e. 
mentioned above, and computed the errors in the following way: 



= ^EAri*sin(0i) 

where R clearly represents the fraction of the error on the transverse direction. 



(14) 



Table 4. Mean values of the magnitude of the position er- 
ror for several tolerances and the transverse component of 
such a vector, for a particular orbit of the perturbed quar- 
tic oscillator at t final = 10^ u.t., written as a percentage. 
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44.3 
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36.6 
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9.352e-4 


46.3 


le-14 


4.929e-8 


31.0 


1.068e-7 


26.3 


9.318e-4 


46.4 


le-15 


3.197e-9 


24.2 


1.068e-7 


26.3 


9.317e-4 


46.4 


le-16 


1.599e-9 


22.2 


1.067e-7 


26.3 


2.543e-2 


50.5 



Table 4 reveals that the major component of the error lies along the trajectory direction. The percent- 
ages for taylor and D0PRI8 can be considered comparable (~ 25 — 30%), although the latter presents 
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less dispersion, while BS gives higher values (~ 45%). 

Taking into account the same orbits as in the previous table, an histogram for each integrator can be 
built, considering e = 10~^^, As shown in Figure [TBI 
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Figure 16. Number of orbits whose fraction of the error in the perpendicular direction to the orbit falls in a given interval 
for taylor, DOPRI8 and BS respectively 



The histograms show that taylor has the majority of its orbits below the mean value, as well as BS, 
while D0PRI8 has a near homogeneity in the interval — 50%. 

In the attempt to determine which of the three integrators is the most precise one, the arithmetic mean 
(for e = 10~^^) of the absolute value of the total error is computed to yield, 

taylor : 3.1971267 x 10"^ 

D0PRI8 : 1.0677529 x 10"^ (15) 
BS : 9.3165436 x 10"^ 

and that of the component in the transverse direction 

{taylor : 7.7557023 x IQ-^O 
D0PRI8 : 2.8064381 x 10"^ (16) 
BS : 4.3187811 x lO""^. 

Therefore, taylor results to be more accurate when computing the coordinates of phase space than 
both D0PRI8 and BS, as already shown in Table 4. In fact, the precision is really good, since it is about 2 
orders of magnitude more precise than D0PRI8 and 6 orders of magnitude more accurate than BS. These 
results support our previous assumption that taylor is a very efficient alternative when dealing with 
problems which involve the solution of the equations of motion and requires the value of the coordinates 
to be quite precise. 
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8. Discussion 

In this work, we reviewed the performances of some variational CIs as well as the numerical techniques for 
their computation over different known dynamical problems. 

As regards the comparison of the variational indicators studied here, we can mention some particular 
points to take into consideration when selecting appropriate tools to study a certain dynamical problem. 
Notwithstanding the excellent performances shown by some techniques, it is advisable to avoid using only 
one indicator of chaos. On the other hand, it is imperative to reduce the number of CIs in the package. 
Therefore, we aimed at selecting a CIsF for a general Hamiltonian system. 

In the previous experiments the RLI and the FLI/OFLI showed the robustest thresholds. The OFLI 
seems to be a reliable suggestion if we can use its time evolution, i.e. to study small samples of orbits or 
for single analysis. The OFLI and the MEGNO are the recommended techniques to identify periodicity 
or stability levels, respecti vely. The GA LI? (or the SALI) -with a proper checking of the i.d.v. because of 
their very low final values, [Barrio . 20051 ]- is also appropria te for identifying periodici t y, and the GALI^. for 



ident if ying regular rnotion on tori of lower dimensionality f jChristodoulidi fc Bouritii . l200()l | . [Skokos et al 



20071 ]. (Skokos e;^ all 12008 ']). Furthermore, to analyze large samples of orbits by means of the final values. 



a combination of the FLI/OFLI and the GALI4 (or the GALl2Ar for a A^-d.o.f. Hamiltonian system) seems 
to be useful to describe the divided phase space of the HH potential. The FLI/OFLI final values did not 
show spurious structures in the regular component within the experiments. The GALI4 with the help 
of the corresponding saturation times is useful both to rough out and describe thoroughly the chaotic 
components. Finally, we ordered the CIs studied according to their CPU-times with a fixed integration 
time. The LI is the most economical in computational time followed closely by the MEGNO and the D is 
the most CPU-time consuming. However, the LI has a rather low speed of convergence. We also showed 
that the FLI/OFLI is less time-consuming due to an efficient implementation of the saturation values. 

The FLI/OFLI proved to be the most versatile technique tested. The FLI showed a robust threshold 
and good performances using the final values to study large samples of orbits and the OFLI showed good 
performances if we use its time evolution for single analysis. Nevertheless, some techniques can be selected 
with the FLI/OFLI to give richer descriptions of the system. 

The previous results showed that a reliable suggestion for a CIsF for the HH potential and for a general 
Hamiltonian fiow is composed of the FLI/OFLI techniques, the MEGNO and the GALl2Ar. Furthermore, 
there is no reason to believe that such CIsF can be only applied to Hamiltonian fiows; it c an also be used 



in mappings, thus improving the efficiency of the CIsF selected in [Maffione et al . 2011a]. Therefore, the 



above-mentioned CIsF can be applied to a general Hamiltonian system. 

We also present the results of a comparative study of the taylor package developed by j Jorba &: Zou 



2OO.5I ] against two other well-known and widely used integrators, namely, D0PRI8 and BS. The three 



methods were applied to diverse dynamical problems and different vector fields as well, with the aim of 
determining their relative efficiency. Both speed and accuracy were tested for three dynamical problems: 
the computation of a fast indicator as the MEGNO, that of the largest LCN (i.e. LI), and the integration 
of the equations of motion for a model of triaxial galaxy and a perturbed quartic oscillator. The numerical 
experiments performed, showed that one of taylor's strongest features is its capacity to preserve the 
energy integral. In this regard, in the more complex system (the triaxial potential), D0PRI8 showed to 
be nearly as accurate as taylor and both turned out to be about two orders of magnitude more accurate 
than BS. In the case of the perturbed quartic oscillator, it resulted about two orders of magnitude more 
precise than both D0PRI8 and BS. 

A characteristic found in this experiments is that the energy preservation yielded by the use of taylor 
turned out to be problem-independent for a given vector field, in contrast to what occured when D0PRI8 
was applied. 

Moreover, the computational effort demanded by the three integrators was subjected to comparison to 
conclude that taylor was considerably slower than D0PRI8 and BS while computing the MEGNO. The 
differences in the required CPU time decreased significantly for the computation of the LI, but even in 
this case taylor remained the slowest one. When dealing with the integration of the equations of motion, 
a change in the behavior was observed, taylor advantaging in speed both D0PRI8 and BS (Figs. [TTland 
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Although this is not a thorough test, in the hght of the results obtained we can say that taylor is 
rather slow when the variational equations or the ones defining the MEGNO are integrated. In general, its 
efficiency seems to diminish as the complexity of the system of ODEs increases. 

However, not only the dynamical problem being attacked, but also the vector field under study impact 
on the differences in the CPU times required by the different integrators. Therefore, for instance, while 
computing the MEGNO the time differences between taylor and the other two integrators are greater for 
the perturbed quartic oscillator than those for the triaxial potential (as seen in Figs. [9land ll2p . 

In an attempt to unravel the observed behavior we adventure to say that whenever dealing with a 
plain formulae for the vector field, such as in the case of the perturbed quartic oscillator, taylor takes too 
much time for computing the recurrences for the automatic differentiation, determining each coefficient and 
building the power series, while the other integrators only evaluate the quite simple functions defining the 
system of ODEs in some points conveniently chosen, thus resulting more efficient than taylor. For a more 
complex system of ODEs, like the ones involved in the model of a triaxial potential, the mere evaluation 
of the field in D0PRI8 or BS would be not as fast as for the quartic oscillator which only encompasses 
polynomic expressions, causing their time differences with taylor to decrease. 

Conversely, regarding the energy preservation the opposite behavior is observed, taylor advantaging 
D0PRI8 for the simpler vector field more than for the one defined by rather cumbersome expressions. In 
fact, in the case of the triaxial potential the obtained differences in /S.Emax between taylor and D0PRI8 
are quite small (barely reaching a factor of ~ 2), while for the perturbed quartic oscillator the errors in 
the energy differ in a factor ranging from ~ 14 to ~ 100 (i.e. in 2 orders of magnitude). Let us say that 
the comparison against BS has been left behind since, even when BS demanded CPU times comparable to 
those required by D0PRI8, the resulting energy preservation in all cases was not as good as that obtained 
by any of the other two integrators. 

In the light of our experiments, a crosswise situation was observed when comparing both the com- 
putational effort and the energy preservation for the integrators being applied to vector fields of different 
complexity. Consequently, generally speaking the supremacy of one integrator over the other could not be 
claimed out. However, what remains clear is that the taylor package results unsuitable for computing the 
MEGNO or integrating the variational equations (this could also be seen on Gerlach & Skokos [2010]), 
where D0PRI8 shows up as a good ch oice to perform the required integrations (let us remark that on the 



3qi 

previous papers [Cincotta et ali l2003l | and [Cincotta &: Simd . l2000l ] . for instance, the D0PRI8 subroutine 



was used for this purpose). 

On integrating only the equations of motion, the analysis of the accuracy reached by means of the 
three integrators exposed that the taylor package turned out to be the most efficient of the three. In this 
application, not only the final errors in the position vector but also their component in the perpendicular 
direction to the orbit were compared in order to estimate how much the trajectory obtained departed from 
the current one. There, taylor was found to be considerably more precise than both D0PRI8 and BS 
(about 2 and 6 orders of magnitude, respectively). 

Thus, considering that the taylor package also revealed itself as advantageous regarding the required 
CPU time and the preservation of the energy integral when integrating the equations of motion of the 3D 
perturbed quartic oscillator, we state that taylor is the most convinient tool to be used in those problems 
that only involve the integration of the equations of motion, such as the study of chaotic diffusion in phase 
space. 

Thanks to investigations like the presented above, we can make smart decisions about the most efficient 
group of CIs and the appropriate numerical integrator to study a given system. Nevertheless, a suitable 
tool should attend the preceeding studies. Thus, in Section [3] we mentioned that the computing of the CIs 
were accomplished with a code we are developing, the LP-VIcode. The LI, the RLI, the MEGNO, the D 
and the SSNs, the SALI and the GALIs and the FLI and the OFLI are already programmed for the systems 
tested, being such systems either a mapping or a flow. They are grouped in several modules, so the user 
can select the desired module to compute the preferred CIs. Moreover, the integration routine used is the 
D0PRI8. However, our final goal is to introduce a code with a growing CIs' database (variational CIs at 
first), where the user can easily program the model to study, and use the desired integrator. 
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